Grizzly bear BC SCR results

library(here)
library(raster)
library(sf)
library(velox)
library(viridis)
library(stringr)
library(rgdal)
library(corrplot)
library(fasterize)
library(ggmap)
library(ggpubr)
library(lubridate)
library(mapview)
library(hrbrthemes)
library(foreach)
library(doParallel)
library(tabularaster)
library(ggrepel)
library(basemaps)
library(RStoolbox)
library(rworldmap)
library(maps)
library(ggsflabel)
library(oSCR)
library(tidyverse)


###set official CRS
off.crs <- "+proj=aea +lat_0=45 +lon_0=-126 +lat_1=50 +lat_2=58.5 +x_0=1000000 +y_0=0 +datum=NAD83 +units=m +no_defs"

##################  
##LOAD DATA
################## 


###Get file paths
filelist = list.files(here::here("Data_Prep","CleanData","CapHists_SECR","CapData"),pattern = ".*.txt", full.names = TRUE)
filelist <- filelist[c(-13,-22, -23, -34,-37:-39)] ###remove sessions with few captures, or primarily rub tree, or outside of extant  grizz range

##load in data
caps <- filelist%>%
  map_df(read_csv)%>%
  as.data.frame()

##spatial data
STACK <- stack(list.files(here::here("Data_Prep", "CleanData", "rasters"),pattern = ".*.tif", full.names = TRUE))
proj4string(STACK) <- off.crs

##rename columns
colnames(caps) <- c("Session", "ID", "Occassion","Detector", "Sex")

##append study area/year onto capture detectors to match trap file
caps <- mutate(caps, Detector=paste(Detector, Session, sep="_"))


######REMOVE RUB TREE TRAP DATA
rtculls<-data.frame()

filelist = list.files(here::here("Data_Prep","CleanData","CapHists_SECR","TrapData"),pattern = ".*.txt", full.names = TRUE)
filelist <- filelist[c(-13,-22, -23, -34,-37:-39)]
filenames <- list.files(here::here("Data_Prep","CleanData","CapHists_SECR","TrapData"),pattern = ".*.txt", full.names = FALSE)
filenames <- filenames[c(-13,-22, -23, -34,-37:-39)]
#filelist <- filelist[c(34)]
trap.list <- list()
n <-c()
for(i in 1:length(filelist)){
  a <- read_csv(filelist[i])
  colnames(a) <- c("Detector", "X", "Y", "Usage_1", "Trap_Type")
  a$sep="/"
  a <- a%>%mutate(Trap_Type=str_split(Trap_Type, "/ ", simplify = TRUE)[,2])%>%
    select("Detector", "X", "Y", "Usage_1", "sep", "Trap_Type")%>%
    mutate(Detector=paste(Detector, str_split(filenames[i], ".txt", simplify = TRUE)[,1], sep="_"))
  
  rtculls <- a%>%filter(Trap_Type%in%"RT")%>%
    select(Detector)%>%
    rbind(rtculls)
  
  a <-  a%>%filter(!Trap_Type%in%" RT")%>%
    select("Detector", "X", "Y", "Usage_1", "sep", "Trap_Type")
  
  use.len <-nchar(a$Usage_1[1])
  a <- tidyr::separate(a,Usage_1, into=paste("Usage",c(1:use.len), sep="_"),sep =c(1:(use.len-1)))
  
  trap.list[[i]]<- as.data.frame(a)
  
  n[i]<-length(unique(a$Detector))
  
  rm(a)
}

caps<- caps%>%
  filter(!Detector %in%rtculls$Detector) ##remove captures at rub trees

###PREP SOME DATA
##occassions per session
k <-caps%>%
  dplyr::group_by(Session)%>%
  dplyr::summarize(sess.max=max(Occassion))%>%
  pull(sess.max)  

##fix up sex column
caps<- caps%>%
  mutate(Sex=case_when(is.na(Sex)~"U", TRUE~Sex))


###SUMMARY STATS
length(unique(caps$ID))
length(unique(caps[caps$Sex%in%"M","ID"]))
length(unique(caps[caps$Sex%in%"F","ID"]))
length(unique(caps[caps$Sex%in%c("U"),"ID"]))
length(unique(caps$Session))
nrow(caps)
length(trap.list)
length(unique(bind_rows(trap.list)$Detector))

##PLOT

##prep data
range <- st_read(here::here("Data_Prep", "Spatial_Layers_ProvSECR", "range.shp"))%>%st_transform(st_crs(STACK))

df <- bind_rows(trap.list)%>%
  select(X,Y, Detector)%>%
  st_as_sf(coords = c("X","Y"),
           crs = off.crs)

##pull city data
data(world.cities)
cities <- world.cities%>%
  filter(lat>46 & lat< 620 & long>(-130) & long<(-110))%>%
  filter(pop>500000|name%in%c("Smithers","Missoula", "Fort Nelson", "Williams Lake"))%>%
  st_as_sf(coords = c("long","lat"),
           crs = 4326)%>%
  st_transform(st_crs(STACK))


##prep basemap
# register_google("xAIzaSyCOwGx2D77XOqRgGhKmcb5F4Kt_S61tCLIx")
# set_defaults(map_service = "osm", map_type = "terrain_bg")
# 
# 
# a <- basemap_raster(st_read("/Users/clayton.lamb/Google Drive/Documents/University/Geographic_Data/Basemaps/bc/shp/province.shp")%>%
#                       st_transform(3857)%>%
#                       st_buffer(400E3))
# a <- a%>%projectRaster(crs=3005)
# writeRaster(a,here::here("Data_Prep","CleanData","basemap","basemap.tif"))
a <- brick(here::here("Data_Prep","CleanData","basemap","basemap.tif"))

world <- getMap(resolution = "high")
world <- st_as_sf(world)
##prep extents
extent = matrix(c(3E5 ,2.5E5 , 3E5 , 17E5 , 20E5  , 17E5, 20E5 ,2.5E5 , 3E5 ,2.5E5  ),ncol=2, byrow=TRUE)
pts = list(extent)
pl1 = st_polygon(pts)
pl1 <- st_sfc(pl1, crs=3005)

inset <- ggplot() +
  geom_sf(data = world,size=0.1, fill="olivedrab", color="black") +
  geom_sf(data=range,alpha=0.5, size=0.2,fill="grey80",inherit.aes = FALSE)+
  geom_sf(data = pl1, fill=NA, color="black", size = 0.7) +
  coord_sf(crs = 3005)+
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.ticks = element_blank(),
        panel.background = element_rect(fill = "slategray3",colour = "black"),
        panel.border = element_rect(fill = NA, color = NA),
        axis.text = element_blank(),
        axis.title = element_blank(),
        plot.margin=unit(c(0,0,0,0),"mm"),
        legend.position = c(0.65,0.075),
        plot.background = element_rect(fill = "transparent",colour = NA))+
  xlim(-1E6,5.3E6)+
  ylim(-2E6,4.5E6)
inset 

map <- ggRGB(a, r=1, g=2, b=3)+
  theme_bw()+
  geom_sf(data=range,alpha=0.1, inherit.aes = FALSE)+
  geom_sf(data=df,size=0.4, pch = 21,alpha=0.6,inherit.aes = FALSE)+
  xlim(3E5,20E5)+
  ylim(2.5E5,17E5)+
  geom_sf(data=cities, inherit.aes = FALSE, color="grey40")+
  geom_sf_text_repel(data=cities, aes(label=name),inherit.aes = FALSE, size=3)+
  annotation_custom(ggplotGrob(inset), xmin =15.5E5, xmax = 21E5, ymin = 12.5E5, ymax = 17.5E5)+
  annotate("text", x = 5.8E5, y = 5E5, label = "Grizzly bear sampling sites", fontface=2)+
  annotate("text", x = 4.3E5, y = 4E5, label =paste( length(unique(caps$Session)) ,"Study-years",sep=" "),size=4)+
  annotate("text", x = 6.7E5, y = 3E5, label =paste(paste(nrow(caps),"Detections from",sep=" "), paste(length(unique(caps$ID)),"Individuals",sep=" "),sep=" "),size=4)+
  theme(axis.title = element_blank())


####map berry 
bec <- st_read(here::here("Data_Prep","Data","berry","bec.shp"))%>%st_transform(st_crs(STACK))

occ <- ggRGB(a, r=1, g=2, b=3)+
  theme_bw()+
  geom_sf(data=bec,size=0.1, pch = 21,alpha=0.3,inherit.aes = FALSE)+
  xlim(3.2E5,18.5E5)+
  ylim(4E5,17E5)+
  geom_sf(data=cities, inherit.aes = FALSE, color="grey40")+
  annotate("text", x = 5.8E5, y = 5E5, label = "Shrub occurrence", fontface=2,size=3)+
  annotate("text", x = 5.6E5, y = 4E5, label = "25,341 plots",size=3)+
  theme(axis.title = element_blank())


prod <- st_read(here::here("Data_Prep","Data","berry","prod.shp"))%>%st_transform(st_crs(STACK))

prod%>%
distinct(geometry)%>%
  tally()

prod.map <- ggRGB(a, r=1, g=2, b=3)+
  theme_bw()+
  geom_sf(data=prod,size=0.5, pch = 21,inherit.aes = FALSE)+
  xlim(3.2E5,18.5E5)+
  ylim(4E5,17E5)+
  geom_sf(data=cities, inherit.aes = FALSE, color="grey40")+
  annotate("text", x = 5.8E5, y = 5E5, label = "Fruit productivity", fontface=2,size=3)+
  annotate("text", x = 5.6E5, y = 4E5, label = "1,566 plots",size=3)+
  theme(axis.title = element_blank())


ggarrange(map,                                                
          ggarrange(occ, prod.map, ncol = 2, labels = c("B", "C")),
          nrow = 2, 
          labels = "A",
          heights=c(2,1))

ggsave(here::here("output", "plots", "map_paper.png"), height=8, width=8, unit="in")
rm(a)
rm(bec)
rm(cities)
rm(inset)
rm(map)
rm(prod)
rm(prod.map)
rm(occ)
rm(world)
rm(world.cities)
rm(extent)
rm(pl1)
rm(pts)
rm(bbox_clip)

######
##add in trap covariates for custom session specs
######

##manually specified
# filelist%>%str_split("/", simplify = TRUE)%>%as_tibble%>%pull(V16) ##get order
# caps%>%group_by(Session)%>%summarize(n=n_distinct(ID,Occassion,Detector))%>%print(n=42) ##see caps, helps to decide groups
d.pool <- factor(c("a1", "a2","a3","a4","a5","a6","a7","a8","a9","a10","Stikine","Stikine", "Nass","Nass","a15","a16","a17","a18","a19","a20",
                   "a21","NorthPurc","SPurc", "a24", "SR","SR","SR","SR","SR","SR","SR","SR","SR","a34","SN","SN","SN","SN","SN","SN","a41","a42","a42","a42"))

ph.pool <- factor(c("Grp1","Grp2","Grp3","Grp4","Grp7","Grp7","SR","KG","Grp1","Stikine",
                    "Stikine","Stikine", "SR","SR","Grp5","Grp5","Grp3","Grp2","PWC","Grp5",
                   "Grp6","Grp1","Grp1", "Grp4", "SR","SR","SR","SR","SR","SR","SR","SR","SR",
                   "PWC","SN","SN","SN","SN","SN","SN","Grp4","Grp2","Grp2","Grp2"))

length(levels(d.pool))

length(levels(ph.pool))

##ADD
for(i in 1:length(trap.list)){
  trap.list[[i]]$d.pool <- d.pool[i]
  trap.list[[i]]$ph.pool <- ph.pool[i]
}



####################################
######GET INTO OSCR FORMAT
####################################
caps$Sex <- as.factor(caps$Sex)
ms.data <- data2oscr(edf=caps, 
                     tdf=trap.list,
                     sess.col=which(colnames(caps) %in% "Session"),
                     id.col=which(colnames(caps) %in% "ID"),
                     occ.col=which(colnames(caps) %in% "Occassion"),
                     trap.col=which(colnames(caps) %in% "Detector"),
                     sex.col=which(colnames(caps) %in% "Sex"),
                     sex.nacode="U",
                     trapcov.names=c("d.pool","ph.pool"),
                     K=k,
                     ntraps=n,
                     tdf.sep="/")



ms.sf <- ms.data$scrFrame

ms.ssDF <- make.ssDF(ms.sf, buffer=35E3, res = 5E3)

## habitat mask
hab <- raster(here::here("Data_Prep", "Spatial_Layers_ProvSECR", "clean", "nonhab.tif"))

###smooth a bit
hab <- focal(hab, w=matrix(1,9,9), fun=mean, na.rm=TRUE,pad=TRUE, padValue=1)%>%
  projectRaster(STACK[[1]])##smooth over ~2500 km2, i.e, the mask size (5x5 km)
hab <-hab>0.6
hab[hab==1] <-NA
hab[hab==0] <-1
plot(hab)

writeRaster(hab,here::here("Data_Prep", "Spatial_Layers_ProvSECR", "hab", "hab.tif"), overwrite=TRUE)

##range mask
range <- st_read(here::here("Data_Prep", "Spatial_Layers_ProvSECR", "range.shp"))%>%
  st_transform(proj4string(STACK[[1]]))%>%
  fasterize(hab, field = "Dis", fun="first")

plot(range)

mask.ms <- (hab*range)
plot(mask.ms)

mask.ms.buff <- hab*(st_read(here::here("Data_Prep", "Spatial_Layers_ProvSECR", "range.shp"))%>%
  st_transform(proj4string(hab))%>%
  st_buffer(15E3)%>%
  fasterize(hab, field = "Dis", fun="first"))

##create function to clip mask of integration by habitat areas only
ms.clipmask <- function(ssDF,raster,p.intact){
  if (!proj4string(raster)%in% c(off.crs, "+proj=aea +lat_0=45 +lon_0=-126 +lat_1=50 +lat_2=58.5 +x_0=1000000 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0")) {
    cat("fix proj to official crs (off.crs)", 
        fill = TRUE)
    return(NULL)
  }
  
  
  
  for(i in 1:length(ssDF)){
    a<-tibble(X=ssDF[[i]][,1],
              Y=ssDF[[i]][,2])%>%
      st_as_sf(coords = c("X","Y"),
               crs = "+proj=aea +lat_1=50 +lat_2=58.5 +lat_0=45 +lon_0=-126 +x_0=1000000 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0")
    
    
    a <- a%>%mutate(cull=velox(raster)$extract_points(a))%>%
      mutate(X=st_coordinates(.)[,1],
             Y=st_coordinates(.)[,2],
             Tr=1)%>%
      as_tibble()%>%
      filter(cull>p.intact)%>%
      select(X,Y,Tr)%>%
      as.data.frame()
    
    
    ssDF[[i]]<-a
  }
  
  return(ssDF)
}

##clip
ms.ssDF <- ms.clipmask(ssDF=ms.ssDF, 
                       raster=mask.ms,
                       p.intact=0.5) ##arbitrary, as the mask is 1/0 so anything in between will split them
#plot(ms.ssDF,ms.sf)

saveRDS(ms.ssDF, here::here("Data_Prep", "Spatial_Layers_ProvSECR","clean","state_space.rds"))

ggplot(data=ms.ssDF[[31]],aes(x=X,y=Y,fill=Tr))+geom_raster()+
  scale_fill_gradientn(colors=rev(viridis_pal()(20)))+theme_minimal()

###cleanup objects
rm(df)
rm(map)
rm(rtculls)
rm(caps)
rm(trap.list)
rm(bbox)
rm(basemap)
rm(ms.data)

###mask extract function
#function
ms.addcovs <- function(ssDF,rasterstack){
  if (!proj4string(rasterstack)%in% c(off.crs, "+proj=aea +lat_0=45 +lon_0=-126 +lat_1=50 +lat_2=58.5 +x_0=1000000 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0")) {
    cat("fix proj to official crs (off.crs)", 
        fill = TRUE)
    return(NULL)
  }
  
  
  for(i in 1:length(ssDF)){
    a<-tibble(X=ssDF[[i]][,1],
              Y=ssDF[[i]][,2])%>%
      st_as_sf(coords = c("X","Y"),
               crs = off.crs)%>%
      as("Spatial")
    
    
    
    ##extract stack
    a <- raster::extract(rasterstack,a, na.rm=TRUE, fun=mean, df=TRUE)
    
    ##add to ssDF
    ssDF[[i]]<-cbind(ssDF[[i]],a%>%select(-ID))
    
    print(paste(i, "of", length(ssDF), "sessions", sep=" "))
    
  }
  
  return(ssDF)
}

#ms.ssDF <- readRDS(here::here("Data_Prep", "Spatial_Layers_ProvSECR","clean","state_space.rds"))

ms.ssDF <- ms.addcovs(ssDF=ms.ssDF,
                      rasterstack=STACK)

saveRDS(ms.ssDF, here::here("Data_Prep", "Spatial_Layers_ProvSECR","clean","state_space_vars.rds"))

ms.ssDF <- readRDS(here::here("Data_Prep", "Spatial_Layers_ProvSECR","clean","state_space_vars.rds"))

##plot
ms.ssDF[[31]]%>%
  select(X, Y, cc_scale, fruit_cal_sum_subset_trim_scale, secure_scale)%>%
  gather("vari", "value", -X,-Y)%>%
  ggplot(aes(x=X,y=Y,fill=value))+geom_raster()+
  scale_fill_gradientn(colors=rev(viridis_pal()(20)))+
  theme_minimal()+
  facet_wrap( ~ vari, ncol=3)

###PLOT
cor.dat <-ms.ssDF
class(cor.dat) <-"list"
M <- cor(dplyr::bind_rows(cor.dat, .id = "column_label")%>%
      select(fruit_cal_sum_subset_trim_scale , vacc_scale, vacc2_scale, f.alt_scale,sprveg_scale, commonveg_scale,ndvi_scale, cc_scale, deltaNDVI_scale , salm_iso_gm_scale,salm_tony_scale, salm_ce_diversity_scale, meat_iso_gm_scale, secure_scale , bb_scale , hum_dens_log_scale, isolation_scale,attractants_log_scale,PINUALB_occcov_scale), use="complete.obs")

corrplot(M, method = "circle", type = "upper", order = "hclust")

##############################
#### POOLED DENSITY SPECIFICATION
##############################
for(i in 1:length(ms.ssDF)){
  ms.ssDF[[i]]$d.pool <- d.pool[i]
    ms.ssDF[[i]]$ph.pool <- ph.pool[i]
}

ms.sf$sigCovs$d.pool <- sapply(ms.ssDF,function(x)x$d.pool[1])
ms.sf$sigCovs$ph.pool <- sapply(ms.ssDF,function(x)x$ph.pool[1])
##summarize captures
ms.sf

################
###start with seeing if it is possible to fit session-specific models, didn't work,
#decided to pool similar sessions
################
sess.pool.mod <- oSCR::oSCR.fit(
        model = list(D ~ d.pool, p0 ~ d.pool + sex + b, sig ~ d.pool + sex),
        scrFrame=ms.sf,
        ssDF=ms.ssDF,
        trimS=35E3)

saveRDS(sess.pool.mod, here::here("output", "mods", "sess.pool.mod.RDS"))

#load models
sess.pool.mod <- readRDS(here::here("output", "mods", "sess.pool.mod.RDS"))

###predict results
#make a dataframe of values predictions
pred.df <- data.frame(session=filelist%>%str_split("/", simplify = TRUE)%>%as_tibble%>%pull(V10)%>%str_sub(1,-5),
                        pool = as.character(d.pool),
                        d.pool = as.character(d.pool),
                      sex=factor(1,levels=c(0,1)),
                      b=0,
                        lab=ph.pool)%>%
  dplyr::mutate(session=case_when(!session%in%c("HWY3_2004","HWY3_2005","Region_2004","Region_2007","S_Purcell_1998","S_Purcell_2001")~str_sub(session,1,-6),
                           TRUE~session%>%as.character()))
#now predict
d.preds <- get.real(model = sess.pool.mod, newdata = pred.df)

sig.preds <- get.real(model = sess.pool.mod, type = "sig", newdata = pred.df)

p0.preds <- get.real(model = sess.pool.mod, type = "det", newdata = pred.df)


##sex column is duplicated, rename
names(d.preds) <- make.unique(names(d.preds))


##drop numeric sex
d.preds <-d.preds%>%select(-sex)%>%rename(sex=sex.1)

#put together
d <-d.preds%>%
  distinct(session,sex, .keep_all = TRUE)%>%
  group_by( session,d.pool,lab)%>%
  summarise_at(vars(estimate,se,lwr,upr), ~ sum(.x)*40)%>%
  left_join(sig.preds%>%select(session, sig=estimate,sig.se=se))%>%
  left_join(p0.preds%>%select(session, p0=estimate,p0.se=se))%>%
  distinct(d.pool, .keep_all=TRUE)

###PLOT with groups
ggplot(d, aes(x=estimate,y=sig))+
  geom_point()+
 theme_ipsum()+
     theme(axis.title.x = element_text(size=15, vjust=-7),
        axis.title.y = element_text(size=15, vjust=2),
        axis.text.x = element_text(size=12),
        axis.text.y = element_text(size=12),
        plot.title = element_text(size=22),
        plot.subtitle = element_text(size=17))+
  labs(title="Density and home range index (Sigma)",
       x="Density", y="Sigma")

ggplot(data=d, aes(x=p0,y=sig, label=session))+
  geom_point(data=d, aes(x=p0,y=sig, label=session, color=lab), size=2)+
  geom_errorbar(data=d, aes(xmin=p0-p0.se, xmax=p0+p0.se, color=lab))+
    geom_errorbar(data=d, aes(ymin=sig-sig.se, ymax=sig+sig.se, color=lab))+
  theme_ipsum()+
  geom_text_repel(segment.size  = 0.2,
    segment.color = "grey50",
    size=3)+
   theme(axis.title.x = element_text(size=15, vjust=-7),
        axis.title.y = element_text(size=15, vjust=2),
        axis.text.x = element_text(size=12),
        axis.text.y = element_text(size=12),
        plot.title = element_text(size=22),
        plot.subtitle = element_text(size=17),
        legend.position = "bottom")+
  labs(title="Sigma vs p0, with detection groups identified",
       x="p0", y="Sigma", color="Detection group")+
  scale_color_brewer(palette="Paired")

lm(sig~estimate, d)%>%summary()

lm(p0~estimate, d)%>%summary()

lm(sig~p0, d)%>%summary()

mods.det <- make.mods(density= c(~d.pool),
                    detection = c(~1,
                                  ~sex,
                                  ~sex + b,
                                  ~ph.pool,
                                  ~ph.pool + sex,
                                  ~ph.pool + sex + b), 
                    sigma = c(~1,
                                  ~sex,
                                  ~ph.pool,
                                  ~ph.pool + sex))

# keep a subset of models to test, needn't be every pairwise option
mods.det <-mods.det[c(1,6,10,15,20,24),]


##run in parallel
cl <- makeCluster(6)    #make the cluster
registerDoParallel(cl)  #register the cluster

mods.det.out <- foreach(i=1:nrow(mods.det),.packages = "oSCR") %dopar% {
  
  m <- list(mods.det[i,1][[1]], # ith model
            mods.det[i,2][[1]],
            mods.det[i,3][[1]],
            mods.det[i,4][[1]]) 
  
  out <- oSCR.fit(m, ms.sf, ssDF=ms.ssDF, trimS=35E3)
  return(out)
}

stopCluster(cl)


 ##examine model fits
(fitList.oSCR(mods.det.out)%>%
  modSel.oSCR())$aic.tab%>%
  as_tibble()

##save outputs
(fitList.oSCR(mods.det.out)%>%
    modSel.oSCR())$aic.tab%>%
  as_tibble()%>%
  write_csv(here::here("output", "tables", "mods.det.mods.csv"))


saveRDS(mods.det.out, here::here("output", "mods", "det.mods.RDS"))

mods.det.out <- readRDS(here::here("output", "mods", "det.mods.RDS"))

###predict results
#make a dataframe of values predictions
pred.df <- data.frame(session=filelist%>%str_split("/", simplify = TRUE)%>%as_tibble%>%pull(V10)%>%str_sub(1,-5),
                        pool = as.character(d.pool),
                        d.pool = as.character(d.pool),
                        ph.pool = as.character(ph.pool),
                      sex=factor(1,levels=c(0,1)),
                      b=0,
                        lab=ph.pool)%>%
  dplyr::mutate(session=case_when(!session%in%c("HWY3_2004","HWY3_2005","Region_2004","Region_2007","S_Purcell_1998","S_Purcell_2001")~str_sub(session,1,-6),
                           TRUE~session%>%as.character()))



#now predict
d.preds <- get.real(model = mods.det.out[[6]], newdata = pred.df)

sig.preds <- get.real(model = mods.det.out[[6]], type = "sig", newdata = pred.df)

p0.preds <- get.real(model = mods.det.out[[6]], type = "det", newdata = pred.df)

##sex column is duplicated, rename
names(d.preds) <- make.unique(names(d.preds))

##drop numeric sex
d.preds <-d.preds%>%select(-sex)%>%rename(sex=sex.1)

#put together
d2 <-d.preds%>%
  distinct(session,sex, .keep_all = TRUE)%>%
  group_by( session,d.pool,lab)%>%
  summarise_at(vars(estimate,se,lwr,upr), ~ sum(.x)*40)%>%
  left_join(sig.preds%>%select(session, sig=estimate,sig.se=se))%>%
  left_join(p0.preds%>%select(session, p0=estimate,p0.se=se))%>%
  distinct(d.pool, .keep_all=TRUE)

###PLOT with groups
ggplot(data=d2, aes(x=p0,y=sig, label=session))+
  geom_point(data=d2, aes(x=p0,y=sig, color=lab), size=2)+
  theme_ipsum()+
  # geom_text_repel(segment.size  = 0.2,
  #   segment.color = "grey50",
  #   size=3)+
   theme(axis.title.x = element_text(size=15, vjust=-7),
        axis.title.y = element_text(size=15, vjust=2),
        axis.text.x = element_text(size=12),
        axis.text.y = element_text(size=12),
        plot.title = element_text(size=22),
        plot.subtitle = element_text(size=17),
        legend.position = "bottom")+
  labs(title="Sigma vs p0, with detection groups identified",
       x="p0", y="Sigma", color="Detection group")+
  scale_color_brewer(palette="Paired")

ggplot(d2, aes(x=estimate,y=sig))+
  geom_point()+
 theme_ipsum()+
     theme(axis.title.x = element_text(size=15, vjust=-7),
        axis.title.y = element_text(size=15, vjust=2),
        axis.text.x = element_text(size=12),
        axis.text.y = element_text(size=12),
        plot.title = element_text(size=22),
        plot.subtitle = element_text(size=17))+
  labs(title="Density and home range index (Sigma)",
       x="Density", y="Sigma")

lm(sig~estimate, d2%>%ungroup%>%distinct(sig,.keep_all = TRUE))%>%summary()

lm(p0~estimate, d2%>%ungroup%>%distinct(sig,.keep_all = TRUE))%>%summary()

lm(sig~p0, d2%>%ungroup%>%distinct(sig,.keep_all = TRUE))%>%summary()

d%>%
  left_join(d2, by="session")%>%
  ggplot(aes(x=estimate.x,y=estimate.y, label=session))+
  geom_point()+
    theme_ipsum()+
  geom_text_repel(segment.size  = 0.2,
    segment.color = "grey50",
    size=3)+
   theme(axis.title.x = element_text(size=15, vjust=-7),
        axis.title.y = element_text(size=15, vjust=2),
        axis.text.x = element_text(size=12),
        axis.text.y = element_text(size=12),
        plot.title = element_text(size=22),
        plot.subtitle = element_text(size=17),
        legend.position = "bottom")+
      geom_abline(intercept = 0, slope = 1)+
  labs(x="Density, detection unpooled", y="Density, detection pooled")

d%>%
  left_join(d2, by="session")%>%
  ggplot(aes(x=se.x,y=se.y, label=session))+
  geom_point()+
    theme_ipsum()+
  geom_text_repel(segment.size  = 0.2,
    segment.color = "grey50",
    size=3)+
   theme(axis.title.x = element_text(size=15, vjust=-7),
        axis.title.y = element_text(size=15, vjust=2),
        axis.text.x = element_text(size=12),
        axis.text.y = element_text(size=12),
        plot.title = element_text(size=22),
        plot.subtitle = element_text(size=17),
        legend.position = "bottom")+
      geom_abline(intercept = 0, slope = 1)+
  labs(x="Density error, detection unpooled", y="Density error, detection pooled")

d%>%
  left_join(d2, by="session")%>%
  ggplot(aes(x=sig.x,y=sig.y, label=session))+
  geom_point()+
    theme_ipsum()+
  geom_text_repel(segment.size  = 0.2,
    segment.color = "grey50",
    size=3)+
   theme(axis.title.x = element_text(size=15, vjust=-7),
        axis.title.y = element_text(size=15, vjust=2),
        axis.text.x = element_text(size=12),
        axis.text.y = element_text(size=12),
        plot.title = element_text(size=22),
        plot.subtitle = element_text(size=17),
        legend.position = "bottom")+
      geom_abline(intercept = 0, slope = 1)+
  labs(x="Sigma, detection unpooled", y="Sigma, detection pooled")

d%>%
  left_join(d2, by="session")%>%
  ggplot(aes(x=p0.x,y=p0.y, label=session))+
  geom_point()+
    theme_ipsum()+
  geom_text_repel(segment.size  = 0.2,
    segment.color = "grey50",
    size=3)+
   theme(axis.title.x = element_text(size=15, vjust=-7),
        axis.title.y = element_text(size=15, vjust=2),
        axis.text.x = element_text(size=12),
        axis.text.y = element_text(size=12),
        plot.title = element_text(size=22),
        plot.subtitle = element_text(size=17),
        legend.position = "bottom")+
      geom_abline(intercept = 0, slope = 1)+
  labs(x="p0, detection unpooled", y="p0, detection pooled")

####WHERE's BURNT
d%>%
  left_join(d2, by="session")%>%
  ggplot(aes(x=sig.se.x,y=sig.se.y, label=session))+
  geom_point()+
    theme_ipsum()+
  geom_text_repel(segment.size  = 0.2,
    segment.color = "grey50",
    size=3)+
   theme(axis.title.x = element_text(size=15, vjust=-7),
        axis.title.y = element_text(size=15, vjust=2),
        axis.text.x = element_text(size=12),
        axis.text.y = element_text(size=12),
        plot.title = element_text(size=22),
        plot.subtitle = element_text(size=17),
        legend.position = "bottom")+
      geom_abline(intercept = 0, slope = 1)+
  labs(x="Sigma error, detection unpooled", y="Sigma error, detection pooled")

d%>%
  left_join(d2, by="session")%>%
  ggplot(aes(x=p0.se.x,y=p0.se.y, label=session))+
  geom_point()+
    theme_ipsum()+
  geom_text_repel(segment.size  = 0.2,
    segment.color = "grey50",
    size=3)+
   theme(axis.title.x = element_text(size=15, vjust=-7),
        axis.title.y = element_text(size=15, vjust=2),
        axis.text.x = element_text(size=12),
        axis.text.y = element_text(size=12),
        plot.title = element_text(size=22),
        plot.subtitle = element_text(size=17),
        legend.position = "bottom")+
      geom_abline(intercept = 0, slope = 1)+
  labs(x="p0 error, detection unpooled", y="p0 error, detection pooled")

mods.d <- make.mods(density= c(~secure_scale + hum_dens_log_scale, ##top down only
                               ~vacc_scale + ndvi_scale + cc_scale + salm_iso_gm_scale + meat_iso_gm_scale, ##bottom up only
                               ~vacc_scale + ndvi_scale + cc_scale + salm_iso_gm_scale + meat_iso_gm_scale + secure_scale + hum_dens_log_scale, ##together
                               ~vacc_scale + ndvi_scale + cc_scale + salm_iso_gm_scale + secure_scale + hum_dens_log_scale, ##together without meat
                               ~vacc_scale + ndvi_scale + cc_scale + secure_scale + hum_dens_log_scale, ##together without salmon
                               ~salm_iso_gm_scale + meat_iso_gm_scale + secure_scale + hum_dens_log_scale, ##meat + top down
                               ~ndvi_scale + cc_scale + secure_scale + hum_dens_log_scale, ##green + top down
                               ~vacc_scale + ndvi_scale + cc_scale + salm_iso_gm_scale + salm_ce_diversity_scale + meat_iso_gm_scale + hum_dens_log_scale, ##together without secure
                               ~vacc_scale + ndvi_scale + cc_scale + salm_iso_gm_scale + secure_scale + hum_dens_log_scale, #together without meat 
                               ~vacc_scale + ndvi_scale + cc_scale + salm_iso_gm_scale + secure_scale + hum_dens_log_scale + bb_scale, #plus bb
                               ~vacc_scale + ndvi_scale + cc_scale + salm_iso_gm_scale + secure_scale + hum_dens_log_scale + isolation_scale, #plus isolation
                               ~vacc_scale + ndvi_scale + cc_scale + salm_iso_gm_scale + secure_scale + hum_dens_log_scale + PINUALB_occcov_scale), #whitebark pine
                    detection = c(~sex + ph.pool + b), 
                    sigma = c(~sex + ph.pool))




##run in parallel
cl <- makeCluster(4)    #make the cluster
registerDoParallel(cl)  #register the cluster

mods.d.out <- foreach(i=1:nrow(mods.d),.packages = "oSCR") %dopar% {
  
  m <- list(mods.d[i,1][[1]], # ith model
            mods.d[i,2][[1]],
            mods.d[i,3][[1]],
            mods.d[i,4][[1]]) 
  
  out <- oSCR.fit(m, ms.sf, ssDF=ms.ssDF, trimS=35E3)
  return(out)
}

stopCluster(cl)


 ##examine model fits
(fitList.oSCR(mods.d.out)%>%
  modSel.oSCR())$aic.tab%>%
  as_tibble()

##save outputs
(fitList.oSCR(mods.d.out)%>%
    modSel.oSCR())$aic.tab%>%
  as_tibble()%>%
  write_csv(here::here("output", "tables", "dens.mods.csv"))


saveRDS(mods.d.out, here::here("output", "mods", "dens.mods.RDS"))

##load
mods.d.out <- readRDS(here::here("output", "mods", "dens.mods.RDS"))
 
##examine model fits
(fitList.oSCR(mods.d.out)%>%
  modSel.oSCR())$aic.tab%>%
  as_tibble()

mods.d.alt <- make.mods(density= c(~fruit_cal_sum_subset_trim_scale + ndvi_scale + cc_scale + salm_iso_gm_scale + secure_scale + hum_dens_log_scale + bb_scale, ##alt fruit
                                   ~f.alt_scale + ndvi_scale + cc_scale + salm_iso_gm_scale + secure_scale + hum_dens_log_scale + bb_scale, ##alt fruit
                                   ~vacc2_scale + ndvi_scale + cc_scale + salm_iso_gm_scale + secure_scale + hum_dens_log_scale + bb_scale, ##alt fruit
                                   
                                   ~vacc_scale + ndvi_scale + cc_scale + salm_f95_Adams_scale + secure_scale + hum_dens_log_scale + bb_scale, ##alt salmon
                                   ~vacc_scale + ndvi_scale + cc_scale + salm_iso_gm_scale + salm_ce_diversity_scale + secure_scale + hum_dens_log_scale + bb_scale, ##plus salmon diversity
                                   ~vacc_scale + ndvi_scale + cc_scale + salm_tony_scale + secure_scale + hum_dens_log_scale + bb_scale, ##alt salmon
                                   
                                   ~vacc_scale + commonveg_scale + cc_scale + salm_iso_gm_scale + secure_scale + hum_dens_log_scale + bb_scale, ##alt green
                                   ~vacc_scale + sprveg_scale + cc_scale + salm_iso_gm_scale + secure_scale + hum_dens_log_scale + bb_scale), ##alt green
                                   
                                   
                                   
                    detection = c(~sex + ph.pool + b), 
                    sigma = c(~sex + ph.pool))




##run in parallel
cl <- makeCluster(4)    #make the cluster
registerDoParallel(cl)  #register the cluster

mods.d.alt.out <- foreach(i=1:nrow(mods.d.alt),.packages = "oSCR") %dopar% {
  
  m <- list(mods.d.alt[i,1][[1]], # ith model
            mods.d.alt[i,2][[1]],
            mods.d.alt[i,3][[1]],
            mods.d.alt[i,4][[1]]) 
  
  out <- oSCR.fit(m, ms.sf, ssDF=ms.ssDF, trimS=35E3)
  return(out)
}

stopCluster(cl)


 ##examine model fits
(fitList.oSCR(mods.d.alt.out)%>%
  modSel.oSCR())$aic.tab%>%
  as_tibble()

##save outputs
(fitList.oSCR(mods.d.alt.out)%>%
    modSel.oSCR())$aic.tab%>%
  as_tibble()%>%
  write_csv(here::here("output", "tables", "dens.alt.mods.csv"))


saveRDS(mods.d.alt.out, here::here("output", "mods", "dens.mods.alt.RDS"))

##load
mods.d.alt.out <- readRDS(here::here("output", "mods", "dens.mods.alt.RDS"))
 
##examine model fits
(fitList.oSCR(mods.d.alt.out)%>%
  modSel.oSCR())$aic.tab%>%
  as_tibble()

##save top mod
top.mod.eco <- mods.d.alt.out[[5]]

##calc linear predictor
predict.denssurf <- function(mod=mod, stack=stack){
  intercept <- mod$coef[mod$coef$param=="d0.(Intercept)","mle"]
  coef <- mod$coef.mle[-c(1:28, nrow(mod$coef.mle)),]  ##need to set 1:however many sig and d0 coeffs there are
  names <- coef[,"param"]
  
  
  coef$match <- c(match(names,paste0("d.beta.",names(stack))))
  
  coef$interact <- str_detect(coef$param,":")%>%as.numeric()
  
  coef$match.x1 <- NA
  coef$match.x2 <- NA
  
  for(i in 1:nrow(coef)){
    if(coef$interact[i]==1){
      names.split <-coef[i,"param"]%>%
        str_remove("d.beta.")%>%
        str_split(":")
      
      
      coef$match.x1[i] <- match(names.split[[1]], names(stack))[1]
      coef$match.x2[i] <- match(names.split[[1]], names(stack))[2]
    }
    
  }
  
  dens.stack <- stack()
  for(i in 1:nrow(coef)){
    
    if(coef$interact[i]==0){
      dens.stack <-addLayer(dens.stack,coef[i,"mle"]*stack[[coef[i,"match"]]])
    }
    
    if(coef$interact[i]==1){
      dens.stack <-addLayer(dens.stack,coef[i,"mle"]*(stack[[coef[i,"match.x1"]]]*stack[[coef[i,"match.x2"]]]))
    }
    
    
  }
  
  dens.stack <-calc(dens.stack, sum)+intercept
  dens.stack <-exp(dens.stack)*0.01
  
  return(dens.stack)
  rm(dens.stack)
}


dens.rast.unclipped <- predict.denssurf(mod=top.mod.eco, stack=STACK)

##aggregate
dens.rast.unclipped <- aggregate(dens.rast.unclipped, 2, fun=sum)


##clip to habitat
range.clip <- mask.ms%>%aggregate(2, fun=max)
dens.rast <- dens.rast.unclipped*range.clip
plot(dens.rast)

sum(values(dens.rast), na.rm=TRUE)

writeRaster(dens.rast, here::here("output", "surface","density_v1.tif"), overwrite=TRUE)

ci <- read.csv(here::here("Data_Prep", "Spatial_Layers_ProvSECR", "CI", "data","CI_MASTER_22Jan2019_CLEAN.csv"))%>%
  select(KILL_CODE, KILL_DATE, SEX, ZONE_NO, GRID_EAST, GRID_NORTH, TOOTH_AGE)%>%
  filter(GRID_EAST>0 & ZONE_NO>0)%>%
  mutate(GRID_EAST=as.numeric(paste0(GRID_EAST,"000")),
         GRID_NORTH=as.numeric(paste0(GRID_NORTH,"000")),
         YEAR=year(ymd(KILL_DATE)),
         count=1)


####Make Bear Data Spatial
Zones<-unique(ci$ZONE_NO)
all <- list()
##Loop to get each zone into correct UTM's and then project to lat long
for(i in 1:length(Zones)){
  a<-ci %>%
    filter(ZONE_NO %in% Zones[i])%>%
    st_as_sf(coords=c("GRID_EAST", "GRID_NORTH"), crs=paste("+proj=utm +zone=",Zones[i]," +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0",sep=""))%>%
    st_transform(crs="+proj=aea +lat_1=50 +lat_2=58.5 +lat_0=45 +lon_0=-126 +x_0=1000000 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0")
  
  
  all[[i]] <- a
}

##Join all CI Zones into a single layer
All<-sf::st_as_sf(data.table::rbindlist(all))
#plot(All)
#mapview(All, zcol="KILL_CODE")



##function to estimate different kill rates (spatial, sexes, and type of kill)
killrate <- function(scale=scale, years=years, sex=sex, killcodes=killcodes, dens=dens){
  res <- res(dens)[1]
  kill.res <- (scale/res)%>%round(0)
  if(kill.res%%2==0){kill.res <- kill.res+1} ##deal with if not odd # (required by raster::focal())
  
  dens.res <- ((scale+6000)/res)%>%round(0) ##add extra to account for home range width
  if(dens.res%%2==0){dens.res <- dens.res+1} ##deal with if not odd # (required by raster::focal())
  
  
  kill <- All%>%
    filter(YEAR%in%years & KILL_CODE%in%killcodes & SEX%in%sex)%>%
    as("Spatial")%>%
    raster::rasterize(y=dens, field = "count", fun=sum)
  
  kill[is.na(kill)] <- 0
  kill <- kill*range.clip
  
  kill <-   (raster::focal(x=kill, w=matrix(1,nrow=kill.res, ncol=kill.res), fun=sum, na.rm=TRUE, pad=TRUE)/length(years))
  
  dens <- dens%>%
    focal(w=matrix(1,nrow=dens.res,ncol=dens.res), fun=sum, na.rm=TRUE, pad=TRUE)
  
  rate <-(kill/dens)*100
  
  if(all(sex==c("F"))){
    rate <-(kill/(dens*0.6))*100
  }
  rate[rate>100] <-100
  
  rate <- projectRaster(rate,STACK[[1]])
  
  return(rate)
}



#Kill Codes:
#1 = hunter,    
#2= animal control
#3= picked-up
#4= Illegal
#5= road mort
#6= rail mort
#9 = trapped

hk.big <- killrate(scale=150E3, years=c(1998:2018), sex=c("F", "M"), killcodes=c(1), dens=dens.rast.unclipped*aggregate(mask.ms.buff,2, fun=max))
nhk.big <- killrate(scale=150E3, years=c(1998:2018), sex=c("F", "M"), killcodes=c(2,4,5,6,9), dens=dens.rast.unclipped*aggregate(mask.ms.buff,2, fun=max))

##make into raster stack
kill.stack <- stack(hk.big, nhk.big)
names(kill.stack) <- c("hk.big", "nhk.big")

plot(kill.stack)

rm(hk.big)
rm(nhk.big)

##scale rasters
layers <- names(kill.stack)
for(i in 1:length(layers)){
  a <- kill.stack[[i]]
  vals <- values(a)[!is.na(values(a))]
  values(a)[!is.na(values(a))] <- (vals-min(vals))/(max(vals)-min(vals))
  names(a) <- paste(layers[i],"scale",sep="_")
  kill.stack <- addLayer(kill.stack,a)
  print(names(a))
}

plot(kill.stack)

##############################
#### SAVE STACK
##############################
writeRaster(kill.stack, filename=here::here("Data_Prep", "CleanData", "rasters.kill",paste0(names(kill.stack),".tif")), bylayer=TRUE,format="GTiff", overwrite=TRUE)

kill.stack <- stack(list.files(here::here("Data_Prep", "CleanData", "rasters.kill"),pattern = ".*.tif", full.names = TRUE))

##############################
##ADD TO STATE SPACE
##############################
ms.ssDF <- readRDS(here::here("Data_Prep", "Spatial_Layers_ProvSECR","clean","state_space_vars.rds"))
ms.ssDF <- ms.addcovs(ssDF=ms.ssDF,
                      rasterstack=projectRaster(kill.stack, crs=off.crs))

ms.ssDF[[31]]%>%
  select(X, Y,paste(c(names(kill.stack)), sep=""))%>%
  gather("vari", "value", -X,-Y)%>%
  ggplot(aes(x=X,y=Y,fill=value))+geom_raster()+
  scale_fill_gradientn(colors=rev(viridis_pal()(20)))+
  theme_minimal()+
  facet_wrap( ~ vari, ncol=4)

##Add rasters to STACK
STACK <- stack(list.files(here::here("Data_Prep", "CleanData", "rasters"),pattern = ".*.tif", full.names = TRUE))
proj4string(STACK) <- off.crs
STACK <- addLayer(STACK, kill.stack)


###PLOT
cor.dat <-ms.ssDF
class(cor.dat) <-"list"
M <- cor(dplyr::bind_rows(cor.dat, .id = "column_label")%>%
      select(fruit_cal_sum_subset_trim_scale , vacc_scale, vacc2_scale, f.alt_scale,sprveg_scale, commonveg_scale,ndvi_scale, cc_scale, deltaNDVI_scale , salm_iso_gm_scale,salm_tony_scale, salm_ce_diversity_scale, meat_iso_gm_scale, secure_scale , bb_scale , hum_dens_log_scale, isolation_scale,attractants_log_scale,PINUALB_occcov_scale, hk.big_scale, nhk.big_scale), use="complete.obs")

corrplot(M, method = "circle", type = "upper", order = "hclust")

##test kill rates
mods.d.kill <- make.mods(density= c(~vacc_scale + ndvi_scale + cc_scale + salm_iso_gm_scale + salm_ce_diversity_scale + secure_scale + hum_dens_log_scale + bb_scale,
                                    ~vacc_scale + ndvi_scale + cc_scale + salm_iso_gm_scale + salm_ce_diversity_scale + secure_scale + hum_dens_log_scale + bb_scale + hk.big_scale,
                                    ~vacc_scale + ndvi_scale + cc_scale + salm_iso_gm_scale + salm_ce_diversity_scale + secure_scale + hum_dens_log_scale + bb_scale + nhk.big_scale,
                                    ~vacc_scale + ndvi_scale + cc_scale + salm_iso_gm_scale + salm_ce_diversity_scale + secure_scale + hum_dens_log_scale + bb_scale + nhk.big_scale + hk.big_scale),
                         detection = c(~sex + ph.pool + b), 
                         sigma = c(~sex + ph.pool))


##run in parallel
cl <- makeCluster(4)    #make the cluster
registerDoParallel(cl)  #register the cluster

mods.kill.out <- foreach(i=1:nrow(mods.d.kill),.packages = "oSCR") %dopar% {
  
  m <- list(mods.d.kill[i,1][[1]], # ith model
            mods.d.kill[i,2][[1]],
            mods.d.kill[i,3][[1]],
            mods.d.kill[i,4][[1]]) 
  
  out <- oSCR.fit(m, ms.sf, ssDF=ms.ssDF, trimS=35E3)
  return(out)
}

stopCluster(cl)


##examine model fits
(fitList.oSCR(mods.kill.out)%>%
    modSel.oSCR())$aic.tab%>%
  as_tibble()


##save outputs
(fitList.oSCR(mods.kill.out)%>%
    modSel.oSCR())$aic.tab%>%
  as_tibble()%>%
  write_csv(here::here("output", "tables", "kill.mods.csv"))


saveRDS(mods.kill.out, here::here("output", "mods", "kill.mods.RDS"))

##load
mods.kill.out <- readRDS(here::here("output", "mods", "kill.mods.RDS"))
 
##examine model fits
(fitList.oSCR(mods.kill.out)%>%
  modSel.oSCR())$aic.tab%>%
  as_tibble()

##top mod
top.mod.pred <- mods.kill.out[[2]]


###predict Density surface
dens.rast2.unclipped <- predict.denssurf(mod=top.mod.pred, stack=STACK)
dens.rast2.unclipped <- aggregate(dens.rast2.unclipped, 2, fun=sum)
dens.rast2 <- dens.rast2.unclipped*range.clip
plot(dens.rast2)

writeRaster(dens.rast2, here::here("output", "surface","density_v2.tif"), overwrite=TRUE)

bc.poly <- st_read("/Users/clayton.lamb/Google Drive/Documents/University/Geographic_Data/Basemaps/bc/shp/province.shp")%>%
  st_transform(proj4string(STACK[[1]]))%>%
  st_simplify(preserveTopology=TRUE, dTolerance=500)

###create BC only mask
bc <- st_read("/Users/clayton.lamb/Google Drive/Documents/University/Geographic_Data/Basemaps/bc/shp/province.shp")%>%
  st_transform(proj4string(STACK[[1]]))%>%
  mutate(Dis=1)%>%
  fasterize(dens.rast2, field = "Dis", fun="first")

plot(bc)

mask <- disaggregate((bc*range.clip),fact=2)
values(mask)[values(mask)<1] <- NA
plot(mask)

##BOTTOM UP

STACK.sens <- subset(STACK, 
                     subset=c("vacc_scale","ndvi_scale","cc_scale","salm_iso_gm_scale","salm_ce_diversity_scale","secure_scale","hum_dens_log_scale","bb_scale"))
names <- names(STACK.sens)
STACK.sens <- STACK.sens*mask
names(STACK.sens)<-names
STACK.sens[["vacc_scale"]] <- STACK.sens[["vacc_scale"]]*1.1
values(STACK.sens[["vacc_scale"]])[values(STACK.sens[["vacc_scale"]])>1]<-1
STACK.sens[["ndvi_scale"]] <- STACK.sens[["ndvi_scale"]]*1.1
values(STACK.sens[["ndvi_scale"]])[values(STACK.sens[["ndvi_scale"]])>1]<-1
STACK.sens[["cc_scale"]] <- STACK.sens[["cc_scale"]]*0.9
STACK.sens[["bb_scale"]] <- STACK.sens[["bb_scale"]]*0.9

dens.snvty.bu <- predict.denssurf(mod=top.mod.eco, stack=STACK.sens)
dens.snvty.bu <- dens.snvty.bu*mask

##TOP DOWN
STACK.sens <- subset(STACK, 
                     subset=c("vacc_scale","ndvi_scale","cc_scale","salm_iso_gm_scale","salm_ce_diversity_scale","secure_scale","hum_dens_log_scale","bb_scale"))
names <- names(STACK.sens)
STACK.sens <- STACK.sens*mask
names(STACK.sens)<-names
STACK.sens[["secure_scale"]] <- STACK.sens[["secure_scale"]]*1.1
values(STACK.sens[["secure_scale"]])[values(STACK.sens[["secure_scale"]])>1]<-1
STACK.sens[["hum_dens_log_scale"]] <- STACK.sens[["hum_dens_log_scale"]]*0.9

dens.snvty.td <- predict.denssurf(mod=top.mod.eco, stack=STACK.sens)
dens.snvty.td <- dens.snvty.td*mask


##BC DENSITY
STACK.sens <- subset(STACK, 
                     subset=c("vacc_scale","ndvi_scale","cc_scale","salm_iso_gm_scale","salm_ce_diversity_scale","secure_scale","hum_dens_log_scale","bb_scale", "hk"))
names <- names(STACK.sens)
STACK.sens <- STACK.sens*mask
names(STACK.sens)<-names
bc.dens <- predict.denssurf(mod=top.mod.eco, stack=STACK.sens)
bc.dens <- bc.dens*mask
sum(values(bc.dens),na.rm=TRUE)

#STACK.sens[["bb_scale"]] <- STACK.sens[["bb_scale"]]*0.99

plot(bc.dens*4000)

##plot changes
plot(dens.snvty.bu-bc.dens)

plot(dens.snvty.td-bc.dens)

##raster math to measure change in density
bu <- ((dens.snvty.bu-bc.dens)/bc.dens)*100
plot(bu)

td <- ((dens.snvty.td-bc.dens)/bc.dens)*100
plot(td)

###PLOT
sens.dat <- rbind(td%>%as.data.frame(xy = TRUE)%>%mutate(Treatment="Top down", estimate=layer)%>%drop_na(estimate)%>%select(x,y,estimate, Treatment),
                  bu%>%as.data.frame(xy = TRUE)%>%mutate(Treatment="Bottom up", estimate=layer)%>%drop_na(estimate)%>%select(x, y,estimate, Treatment))

source("https://gist.githubusercontent.com/benmarwick/2a1bb0133ff568cbe28d/raw/fb53bd97121f7f9ce947837ef1a4c65a73bffb3f/geom_flat_violin.R")

lb <- function(x) mean(x) - sd(x)
ub <- function(x) mean(x) + sd(x)

sumld <- sens.dat %>% 
  select(Treatment,estimate)%>%
  group_by(Treatment) %>%
  summarise_all(funs(mean, median, lower = lb, upper = ub))

sumld

sens.plot <-ggplot(data=sens.dat, aes(x=Treatment, y=estimate, fill=Treatment))+
  #geom_flat_violin(position = position_nudge(x = .2, y = 0), alpha = .8) +
  geom_point(data=sens.dat%>%group_by(Treatment)%>%sample_n(1000),aes(color=Treatment),position = position_jitter(width = .15), size = .3, alpha = 0.1)+
  geom_boxplot(width = .1, outlier.shape = NA, alpha = 0.5) +
  xlab("Treatment")+
  ylab((bquote(paste("Change in bear density (%)"))))+
  scale_fill_viridis_d()+
  scale_color_viridis_d()+
  theme_ipsum()+
  theme(axis.title.x = element_text(size=15, vjust=-7),
        axis.title.y = element_text(size=15, vjust=2),
        axis.text.x = element_text(size=12),
        axis.text.y = element_text(size=12),
        plot.title = element_text(size=22),
        plot.subtitle = element_text(size=17),
        legend.position = "none")

sens.plot

#ggsave(here::here("output", "plots", "limiting_influence.png"), height=4, width=4, unit="in")


###MAP
plot((bu-td)>0)

plot((bu-td)<0)

bu.plot.sens <- sens.dat%>%
  filter(Treatment%in%"Bottom up")%>%
  ggplot(data = .)+
  geom_tile(aes(x = x, y = y, fill =estimate))+
  geom_sf(data=bc.poly, fill=NA,lwd=0.2)+
  scale_fill_viridis(direction=-1)+
  theme_bw()+
  labs(title="Bottom up", fill="%")+
  theme(plot.title = element_text(face = "bold",
                                  size = rel(1.2), hjust = 0.5),
        text = element_text(),
        panel.background = element_rect(colour = NA),
        plot.background = element_rect(colour = NA),
        panel.border = element_rect(colour = NA),
        axis.title=element_blank(),
        axis.text.y = element_text(size = rel(1.1)), 
        axis.text.x = element_text(size = rel(1.1), angle=45, hjust=1),
        axis.line = element_line(colour="black"),
        axis.ticks = element_line(),
        legend.position = c(0.9,0.7),
        panel.grid.major = element_line(colour="#f0f0f0"),
        panel.grid.minor = element_blank(),
        legend.key = element_rect(colour = NA),
        strip.background=element_rect(colour="#f0f0f0",fill="#f0f0f0"),
        strip.text = element_text(face="bold"))



td.plot.sens <- sens.dat%>%
  filter(Treatment%in%"Top down")%>%
  ggplot(data = .)+
  geom_tile(aes(x = x, y = y, fill =estimate))+
  geom_sf(data=bc.poly, fill=NA,lwd=0.2)+
  scale_fill_viridis(direction=-1)+
  theme_bw()+
  labs(title="Top down", fill="%")+
  theme(plot.title = element_text(face = "bold",
                                  size = rel(1.2), hjust = 0.5),
        text = element_text(),
        panel.background = element_rect(colour = NA),
        plot.background = element_rect(colour = NA),
        panel.border = element_rect(colour = NA),
        axis.title=element_blank(),
        axis.text.y = element_text(size = rel(1.1)), 
        axis.text.x = element_text(size = rel(1.1), angle=45, hjust=1),
        axis.line = element_line(colour="black"),
        axis.ticks = element_line(),
        panel.grid.major = element_line(colour="#f0f0f0"),
        panel.grid.minor = element_blank(),
        legend.key = element_rect(colour = NA),
        legend.position = c(0.9,0.7),
        strip.background=element_rect(colour="#f0f0f0",fill="#f0f0f0"),
        strip.text = element_text(face="bold"))



ggarrange(sens.plot, bu.plot.sens,td.plot.sens,
          nrow=1, ncol=3,
          labels="AUTO",
          widths=c(1,1.2,1.2))

ggsave(here::here("output", "plots", "elas.plot_plusmaps.png"), height=4, width=11, unit="in")

coefs <-c("vacc_scale","ndvi_scale","cc_scale","salm_iso_gm_scale","salm_ce_diversity_scale","secure_scale","hum_dens_log_scale","bb_scale")
coefs_raw <-c("vacc","ndvi","cc","salm_iso_gm","salm_ce_diversity","secure","hum_dens_log", "bb")

STACK.pred <- subset(STACK, 
                     subset=coefs)%>%
  tabularaster::as_tibble(cell=FALSE,xy=TRUE)%>%
  drop_na(cellvalue)%>%
  left_join(tibble(dimindex=1:length(coefs), var=coefs))

var.mean <- STACK.pred%>%
  group_by(var)%>%
  summarise(mean=median(cellvalue))



STACK.pred  <- STACK.pred%>%
  mutate(round=round(cellvalue,1)%>%as.character())%>%
  group_by(var,round)%>%
  sample_n(size=3,replace=FALSE)%>%
  left_join(subset(STACK, 
                   subset=coefs_raw)%>%
              tabularaster::as_tibble(cell=FALSE,xy=TRUE)%>%
              drop_na(cellvalue)%>%
              select(raw=cellvalue, x,y,dimindex),
            by=c("dimindex","x","y"))


STACK.pred <- STACK.pred%>%
  ungroup()%>%
  select(var,cellvalue,x,y)%>%
  pivot_wider(values_from=cellvalue, names_from=var)%>%
  left_join(STACK.pred%>%
          ungroup()%>%
          select(var,raw,x,y)%>%
          mutate(id=rep(1:(11*3),times=length(coefs)))%>%
          mutate(var=str_split(var,"_scale", simplify = TRUE)[,1]),
          by=c("x","y"))

##fill NA's
STACK.pred <- STACK.pred%>%
  mutate(hum_dens_log_scale=case_when(is.na(hum_dens_log_scale)~var.mean%>%filter(var%in%"hum_dens_log_scale")%>%pull(mean),
                              TRUE~hum_dens_log_scale),
         ndvi_scale=case_when(is.na(ndvi_scale)~var.mean%>%filter(var%in%"ndvi_scale")%>%pull(mean),
                              TRUE~ndvi_scale),
         salm_ce_diversity_scale=case_when(is.na(salm_ce_diversity_scale)~var.mean%>%filter(var%in%"salm_ce_diversity_scale")%>%pull(mean),
                              TRUE~salm_ce_diversity_scale),
         salm_iso_gm_scale=case_when(is.na(salm_iso_gm_scale)~var.mean%>%filter(var%in%"salm_iso_gm_scale")%>%pull(mean),
                              TRUE~salm_iso_gm_scale),
         secure_scale=case_when(is.na(secure_scale)~var.mean%>%filter(var%in%"secure_scale")%>%pull(mean),
                              TRUE~secure_scale),
         vacc_scale=case_when(is.na(vacc_scale)~var.mean%>%filter(var%in%"vacc_scale")%>%pull(mean),
                              TRUE~vacc_scale),
         cc_scale=case_when(is.na(cc_scale)~var.mean%>%filter(var%in%"cc_scale")%>%pull(mean),
                              TRUE~cc_scale),
         bb_scale=case_when(is.na(bb_scale)~var.mean%>%filter(var%in%"bb_scale")%>%pull(mean),
                              TRUE~bb_scale)
         )
  



pred.dat <- get.real(top.mod.eco, type = "dens", newdata = data.frame(STACK.pred), d.factor = 40)


pred.dat%>%
  mutate(raw=case_when(var%in%"hum_dens_log"~exp(raw)-100,
                       TRUE~raw))%>%
  mutate(cull=case_when(var%in%"vacc" & raw>6500~1,
                       TRUE~0))%>%
  filter(cull==0)%>%
  distinct(var,raw,.keep_all = TRUE)%>%
  group_by(var,raw)%>%
  summarise_at(c("estimate", "se"), sum)%>%
  left_join(tibble(var=coefs_raw, name=c("Berries","Greenness","Canopy cover","Salmon diet","Salmon diversity","Secure habitat","Human density", "Black bear")))%>%
  ggplot(aes(x=raw, y=estimate))+
  geom_line()+
  geom_ribbon(aes(ymin = estimate-se, ymax = estimate+se), fill = "grey70", alpha=0.5)+
  facet_wrap(~name, scales="free")+
  xlab("")+
  ylab((bquote(paste("Bear density (/1000 km"^{2},")"))))+
  scale_fill_viridis_d()+
  theme_bw()+
  theme(plot.title = element_text(face = "bold",
                                              size = rel(1.2), hjust = 0.5),
                    text = element_text(),
                    panel.background = element_rect(colour = NA),
                    plot.background = element_rect(colour = NA),
                    panel.border = element_rect(colour = NA),
                    axis.title = element_text(size = rel(1.3)),
                    axis.title.y = element_text(angle=90,vjust =2),
                    axis.title.x = element_text(vjust = -0.2),
                    axis.text = element_text(size = rel(1.1)), 
                    axis.line = element_line(colour="black"),
                    axis.ticks = element_line(),
                    panel.grid.major = element_line(colour="#f0f0f0"),
                    panel.grid.minor = element_blank(),
                    legend.key = element_rect(colour = NA),
                    legend.position = "none",
                    plot.margin=unit(c(10,5,5,5),"mm"),
                    strip.background=element_rect(colour="#f0f0f0",fill="#f0f0f0"),
                    strip.text = element_text(face="bold"))

ggsave(here::here("output", "plots", "individ_response.png"), height=5, width=6, unit="in")

hab5K <- projectRaster(hab,STACK[[1]])
hab5K <- hab5K%>%aggregate(fact=10, fun=mean, na.rm=TRUE)
hab5K <- hab5K >0.5
hab5K[values(hab5K)<1] <- -1
names(hab5K) <- "hab5K"
coefs <-c("vacc_scale","ndvi_scale","cc_scale","salm_iso_gm_scale","salm_ce_diversity_scale","secure_scale","hum_dens_log_scale","bb_scale","hk.big_scale")
STACK.pred  <- (subset(STACK, 
       subset=coefs))%>%
  raster::aggregate(fact=10, fun=mean, na.rm=TRUE)%>%
  stack(hab5K)%>%
  tabularaster::as_tibble(cell=FALSE,xy=TRUE)%>%
  drop_na(cellvalue)%>%
  left_join(tibble(dimindex=1:(length(coefs)+1), var=c(coefs,"hab5K")))%>%
  select(x,y,var,cellvalue)%>%
  pivot_wider(values_from=cellvalue, names_from=var)%>%
  filter(hab5K>=0)%>%
  mutate(#nhk.big_scale=case_when(is.na(nhk.big_scale)~median(nhk.big_scale,na.rm=TRUE), TRUE~nhk.big_scale),
         hk.big_scale=case_when(is.na(hk.big_scale)~median(hk.big_scale,na.rm=TRUE), TRUE~hk.big_scale))%>%
  drop_na()


##run in parallel
cores = 4 

STACK.pred.split <- STACK.pred %>% 
  group_by((row_number()-1) %/% (n()/cores)) %>%
  nest%>%pull(data)

cl <- makeCluster(cores)    #make the cluster
registerDoParallel(cl)  #register the cluster

bc.dat <- foreach(i=1:cores,.packages = "oSCR") %dopar% {
  dat <-as.data.frame(STACK.pred.split[[i]])
  pred <- get.real(top.mod.eco, type = "dens", newdata = dat)
  #pred <- get.real(top.test, type = "dens", newdata = dat)
  #pred <- get.real(top.mod.eco.nob, type = "dens", newdata = dat)
  return(pred)
}

stopCluster(cl)

bc.dens.mf <- bc.dat%>%
  bind_rows

nrow(bc.dens.mf)/2